Thermal destruction of chiral order in a two-dimensional model of coupled trihedra 
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We introduce a minimal model describing the physics of classical two-dimensional (2D) frustrated 
Heisenberg systems, where spins order in a non-planar way at T — 0. This model, consisting of 
coupled trihedra (or Ising-RP 3 model), encompasses Ising (chiral) degrees of freedom, spin- wave 
excitations and Z2 vortices. Extensive Monte Carlo simulations show that the T — chiral order 
disappears at finite temperature in a continuous phase transition in the 2D Ising universality class, 
despite misleading intermediate-size effects observed at the transition. The analysis of configura- 
tions reveals that short-range spin fluctuations and Z2 vortices proliferate near the chiral domain 
walls explaining the strong renormalization of the transition temperature. Chiral domain walls can 
themselves carry an unlocalized Z2 topological charge, and vortices are then preferentially paired 
with charged walls. Further, we conjecture that the anomalous size-effects suggest the proximity of 
the present model to a tricritical point. A body of results is presented, that all support this claim: 
(i) First-order transitions obtained by Monte Carlo simulations on several related models (ii) Ap- 
proximate mapping between the Ising-RP 3 model and a dilute Ising model (exhibiting a tricritical 
point) and, finally, (iii) Mean-field results obtained for Ising-multispin Hamiltonians, derived from 
the high-temperature expansion for the vector spins of the Ising-RP 3 model. 

PACS numbers: 05.50. +q,75.10.Hk 



I. INTRODUCTION 

On bipartite lattices, the energy of the classical Heisen- 
berg, as well as XY, antiferromagnet is minimized by 
collinear spin configurations. Any two such ground states 
can be continuously transformed into one another by a 
global spin rotation. By contrast, it is quite common 
that the ground state manifold of frustrated magnets 
comprises several connected components, with respect to 
global spin rotations, that transform into one another un- 
der discrete symmetry only. Examples include Villain's 
fully frustrated XY mgdeljil the J\ — J2 Heisenberg model 
on the square, lattice,!™ the J — K4 model on the triang 
gular lattice,cl the J\ — J3 model on the squase, lattice,tl 
and the J\ — J2 model on the kaaome latticeHD 

The Mermin- Wagner theoremta forbids the sponta- 
neous breakdown of continuous symmetries, such as spin 
rotations, at any T > in two dimensions (2D). How- 
ever, as was first noticed by Villainjil the breakdown of 
the discrete symmetries relating the different connected 
components of the ground state manifold may indeed 
give rise to finite temperature phase transition(s). Such 
transitions have been evidenced numerically in a num- 
ber oj^fgugJjrated systems, with either XYIj or Heisenberg 

We are interested in a particular class of models with 
Heisenberg spins ^where the ground state has non-planar 
long-range order .dDE In this case the ground state is la- 
beled by an 0(3) matrix. Hence the ground state man- 
ifold is 0(3) = SO(3) x 7L,2 which breaks down into 
two copies of SO (3) . The two connected components, 
say 1 and 2, are exchanged by a global spin inversion 



(Si — ► —Si) and may be labeled by opposite scalar chi- 
ralities Si ■ (Sj x Sk)- Hence we introduce a local Ising 
variable a(f) = ±1 which measures whether the spins 
around f have the chirality of sector 1 or 2. At T = the 
chiralities a(r) are long-range ordered and the ground 
state belongs to a given sector. On the other hand, at 
high enough temperature the system is fully disordered. 
Hence, on very general grounds we expect the sponta- 
neous breakdown of the spin inversion symmetry, associ- 
ated to (er(r)) 7^ 0, at some intermediate temperature. 

Further, from the standpoint of Landau-Ginzburg the- 
ory, one anticipates a critical transition in the 2D Ising 
universality class- However, of the two relevant models 
studied so far,crQiD none shows the signature of an Ising 
transition. Instead, as was pointed out by some of us in 
Ref. 0, the existence of underlying (continuous) spin de- 
grees of freedom complicates the naive Ising scenario, and 
actually drives the chiral transition towards first-order. 

To get a better sense of this interplay between discrete 
and continuous degrees of freedom, it is useful to remem- 
ber that in 2D, although spin- waves disorder the spins 
at any T > 0, the spin-spin correlation length may be 
huge (£ ~ cxp(— A/T)) at low temperature,t3 especially 
in frustrated systems. Hence, it is likely that the effec- 
tive, spin-wave mediated, interaction between the emer- 
gent Ising degrees of freedom extends significantly be- 
yond one lattice spacing, even at finite temperature and 
in 2D. 

Further, we point out that the excitations built on the 
continuous degrees of freedom are not necessarily lim- 
ited to spin-waves. To be more specific, if the connected 
components of the ground state manifold are not simply 
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connected, as is the case for 50(3), then there also ex- 
ists defects in the spin textures. Here, IIi(50(3)) = Z 2 
implies that Z 2 point defects (vortices in 2D) are topo- 
logical^ stable. Clearly these additional excitations may 
also affect the nature of the transition associated to the 
Ising degrees of freedom. In fact, it was shown on one 
exampleO that the first-order chiral transition is triggered 
by the proliferation of these defects. 

In this paper we aim at clarifying the nature of the in- 
terplay between the different types of excitations found 
in Heisenberg systems with non-planar long-range order 
at T = 0. Note that the associated unit cell is typi- 
cally quite large, which severely limits the sample sizes 
amenable to simulations. Hence, we introduce a minimal 
model with the same physical content as the frustrated 
models studied in Refs ^§0. 

As was already mentioned, in the above frustrated spin 
systems, the spin configuration at T = is entirely de- 
scribed by an 0(3) matrix, or equivalently a trihedron 
in spin space. At low temperature, the spin long-range 
order is wiped out by long wavelength spin waves, but 
from the considerations above we anticipate that, at low 
enough temperatures, the description in terms of tri- 
hedra in spin space still makes sense, at least locally. 
To be more specific, we assign three unit vectors Si 
(a = 1,2,3) to every site i of the square lattice, subject 
the orthogonality constraint 



the interplay between discrete and continuous degrees of 
freedom we perform a microscopic analysis of typical con- 
figurations. 

For the most part, the remainder of our work originates 
from the observation of peculiar intermediate size effects 
at the transition. This leads us to argue that the present 
model lies close to a tricritical point in some parameter 
space. 

To support our claim we first introduce and perform 
Monte Carlo simulations on two modified versions of 
our model that i) preserve the 0(3) manifold of ground 
states, and ii) lead to a first-order transition of the Ising 
variables. 



This is further elaborated on in section |V|, where we 
draw an analogy between the Ising-RP 3 model and the 
large q Potts model. Another analogy, this time to a 
dilute Ising model, is drawn in section |V|, where we argue 
that the regions of strong misalignment of the trihedra, 
near Ising domain walls, can be treated as "depletions" 
in the texture formed by the 4D-vectors. 

Finally, in Section [V| we take another route and trace 
out the continuous degrees of freedom perturbatively, re- 
sulting in an effective model for the Ising variables, that 
we proceed to study at the mean-field level. 



II. BASICS OF THE MODEL 
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and we assume the following interaction energy 
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a=l 



(1) 
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Hence we consider three ferromagnetic Heisenberg mod- 
els, tightly coupled through the rigid constraint (|l|). At 
T = 0, the energy is minimized by any configuration with 
all trihedra aligned, and the manifold of ground states 
is 0(3), as desired. This alone ensures the existence of 
the three types of excitations: i) SO(3) spin waves (cor- 
responding to the rotation of the trihedra), ii) 5*0(3) 
vortices, and hi) Ising (chiral) degrees of freedom, comes 
sponding to the right or left-handedness of the trihedral 
The present study is devoted to the model defined by 
Eqs. | and |. 

In Section [0] we reformulate this model more conve- 
niently in terms of chiralities and four-dimensional (4D) 
vectors, yielding the so-called Ising-RP 3 model. For clar- 
ity we first consider a simplified version of this model 
where the chirality variables are frozen, and we use this 
setup to detail our method to detect vortex cores and 
analyze their spatial distribution. 

In section III, we return to the full Ising-RP 3 model 
and evidence the order-disorder transition of the Ising 
variables at finite-temperature. The nature of the tran- 
sition is asserted by a thorough finite-size analysis us- 
ing Monte-Carlo simulations. To clarify the nature of 



A. Ising-RP formulation 

The model defined by Eqs. [j] and || can be conve- 
niently reformulated as an Ising model coupled to a four- 
component spin system with biquadratic interactions. In- 
deed, every trihedron is represented by an SO(3) matrix 
Mi and a chirality Oi = ±1: 
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Using these variables, Eq. || reads 

> ; a^-Tr [{MiYMj] 



(3) 



{i,0) 



The isomorphism between SO(i) and SU(2)/ {1, — 1}, 
maps a rotation M(9, n) of angle 9 about the n axis 
onto the pair of SU(2) matrices ±exp(i|n ■ a), where 
the components of a are the Pauli matrices. The lat- 
ter can be written using two opposite 4D real vectors 
±v = ±(vo,v x ,v y ,v z ) with v 2 = 1: 



exp [ i—n ■ a 



vq + iv z iv x + v v 



v y v 



IV, 



v = cos(0/2) , v a = n o sin(0/2). 



(5) 
(6) 



Irrespective of the local (arbitrary) choice of represen- 



tation ±Vi, one has 



Tr [(MtfMj] =±{v i -v j y 



(7) 
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FIG. 1: (Color online). SO(3) can be represented by a ball 
of radius n. A rotation of angle 9 E [0, n] around the n axis 
is represented in this ball by the extremity of the vector On. 
Two opposite points of the surface (joined by a dashed line) 
represent the same rotation. A loop with an odd number of 
crossings, such as ABFEA or CDEFC, cannot be continu- 
ously shrunk to a point. This is not the case when the number 
of crossings is even: a) ABCDEA crosses the surface both at 
(AB) and (CD). b)-c) It can be continuously deformed to col- 
lapse points B on C, and A on D. d) AD and BC become 
contractible closed loops, that collapse on identical points. 



whence the energy reads 

E = - J2 Wi ( 4 (^ ' vj? - 1) • (8) 



B. Z2 vortices in the fixed-chirality limit 

Here we first consider the simplified case where the 
Ising degrees of freedom are frozen to, say, <r, =f-i+l- In 
this limit we recover the so-called MP 3 modeljla which 
contains both spin- waves and Z2 vortices, and describes 
interacting 50(3) matrices. Together with related frus- 
trated spin models (Heisenberg model on the triangular 
lattice for instance) , the KLP 3 model has been the central 
subject of a number of studies focusing on a putative 
binding-unbiBd^ff-taBsitiott^|4h^|Zo— vortices at finite 
temperatureJliliMEMliyM.E&^ This is a deli- 
cate and controversial issue, which is not essential to the 
present work. Instead, we merely discuss some proper- 
ties of the vortex configurations that will be useful for 
comparison with the full 0(3) (or Ising-IRP 3 ) model. 

To locate the topological point defects (vortex 
cores) we resort to the usual procedure: Consider 
a closed path, C, on the lattice, running through 
sites io, ■ ■ ■ , i n — i) in = io- £ induces a loop Cc in 
the order parameter space, defined by the matrices 
Mj , ■ • ■ , Mi n _ ± , Mi n (in this space the path between Mi k 
and M ik+1 is defined as the "shortest one"). The homo- 
topy class of Cc is an element of IIi . If it is the identity, 
then the Cc loop is contractible. .jQtherwise, L surrounds 
(at least) one topological defect.c3 

In the RP 3 model, the fundamental group is 
111 (50(3)) = Z2, so that the topological charge can take 
only two values: the homotopy class of Cc correspond 
to the parity of the number of point defects enclosed in 
C This is to be contrasted with the better-known £0(2) 
vortices, associated, for instance, to the XY model: the 



latter carry an integer charge (Hi (50(2)) = Z) while the 
former carry merely a sign. 

The number of Z 2 vortices of a given configuration 
is obtained by looking for vortex cores on each elemen- 
tary plaquette of the lattice. The vorticity f2(p) of a 
square plaquette p is computed using by mapping 50(3) 
to MP 3 = 5 3 /Z2: on every site i we arbitrarily choose 
one of the two equivalent representations ±Vi of the local 
50(3) matrix and compute 

tt{p) = IJsign^.^), (9) 

□ P 

where i and j are nearest neighbors and the product runs 
over the four edges of D p . The associated closed loop 
in 50(3) is non-contractible when the plaquette hosts a 
vortex core, and is identified by 0(p) = — 1. Note that 
is a "gauge invariant" quantity, i.e. it is independent 
of the local choice of representation ±Vi, as it should be. 

We performed a Monte Carlo simulation of the MP 3 
model using a Wang-Landau algorithm, detailed in Ap- 
pendix |A[ In the upper panel of Fig. ^ we plot the vortex 
density nn, defined as the number of plaquettes host- 
ing a vortex core divided by the total number of plaque- 
ttes. Vortices are seen to appear and the density increases 
upon increasing the temperature from T ~ 1. However, 
no critical behavior (scaling) is observed upon increasing 
the system size. The latter is also true of other sim- 
ple thermodynamic quantities, such as the energy or the 
specific heat, although the latter is maximum when the 
increase in nn is steepest, at T ~ 1.3.E3 

A typical configuration is shown in the lower panel of 
Fig. at a temperature T = 1.1, about 20% lower than 
that of the maximum of the specific heat. Plaquettes 
hosting a vortex core are indicated by a (red) bullet. Note 
that nn is rather small at this temperature, and that all 
vortices are paired, except for two, evidencing the strong 
binding of the vortices. 

In the same figure the width of every bond is pro- 
portional to Bij = 1 — (ili ■ vj) 2 and indicates the relative 
orientation of the two 4D-vectors Vi and Vj . B^ = (no 
segment) corresponds to parallel 4D-vectors (or identical 
trihedra), which minimizes the bond energy ( E- t j = — 3). 
In particular, all bonds have Bij = at T — 0. On 
the contrary, Bij = 1 (thick black segment) indicates a 
maximally frustrated bond (Eij = +1) with orthogonal 
4D- vectors (the two trihedra differ by a rotation of angle 
7r) . Figure || shows that the vortex cores are located in 
regions of enhanced short-range fluctuations of the con- 
tinuous variables (represented by thick bonds), but that 
the converse is not necessarily true. 



III. NUMERICAL SIMULATIONS OF THE 
ISING-RP 3 MODEL 

We now return to the full Ising-IRP 3 model, defined 
in Eq. ||. We report results of Monte Carlo simulations, 
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FIG. 2: (Color onlin e) Uniform frozen chirality model (see 
text paragraph [il B| ). Top: Vortex density nn versus tem- 
perature T, from L = 16 to L = 64, from bottom to top. 
Bottom: Thermalized configuration at T ~ 1.1 on a sys- 
tem of size L = 72 (only a 45x45 square snippet is shown). 
Square plaquettes carrying a 1i vortex core are denoted by a 
(red) bullet. The width of each bond is proportional to 
Bij — 1 — (vi ■ Vj) 2 (by slices of 1/8). Bonds with By < 1/8 
are not drawn. 



using the Wang-Landau algorithm described in Appendix 
[a|, for linear sizes up to L = 88 with periodic boundary 
conditions. 



A. Specific heat and energy distribution 

Once the Ising degrees of freedom are relaxed, the max- 
imum of the specific heat diverges with the system size, 
indicating a phase transition. Further, the scaling as 
~ log(L) (Fig. ||) for the largest L suggests a continuous 
transition in the Ising-2D universality class. To ascer- 
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FIG. 3: (Color online) Maximum value of the specific heat as a 
function of system size L. The scaling at large L is correctly 
reproduced by the affine log fit A + Bln(L) (straight line). 
Inset: specific heat C v versus temperature T for L = 64 and 
80 from right to left. 
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FIG. 4: (Color online) Probability distribution P(E, T c ) of the 
energy per site at T C (L) for increasing linear L, from L = 16 
to L = 88, from bottom to top. T C (L) is the temperature 
where C v is maximum. A double-peak is visible for L < 60, 
and disappears for L = 72 and L = 80. 



tain the continuous nature of the transition, we computed 
the probability distribution P(E, T = T c ) of the energy 
per site, obtained from the density of states g(E) (Fig. Q). 
Here T C (L) is defined as the temperature where the spe- 
cific heat is maximum. Rather surprisingly, from small 
to intermediate lattice sizes, P(E,T = T c ) shows the bi- 
modal structure characteristic of a first-order transition. 
However, this feature disappears smoothly upon increas- 
ing the system size, and a single peak finally emerges for 
L > 72. 

Hence we claim that the phase transition is continu- 
ous indeed, in the Ising-2D universality class, although 
the scaling and energy distribution at moderate sizes is 
misleading. This peculiar finite-size behavior evidences a 
large but finite length scale whose exact nature has not 
been elucidated so far. We conjecture that it may be 
related to the proximity, in some parameter space, to a 
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tricritical point where the transition becomes discontin- 
uous. Further argum ents in support to this claim will be 
provided in Sections HIE , IV and 0. 



B. Ising order parameter 

The ordering of the chirality variables <7j is probed by 
the following Ising order parameter: 



(10) 



where N is the number of sites. With the results above 
in mind we analyze the finite-size effects on the chirality 

using the scaling law L^ u a = f [l x I v (t^)) with 

(3 = 1/8 and v = 1. Fig. 5(a) shows a data collapse of 
the Ising order parameter for 40 < L < 88, in agreement 
with th e 2D -Ising critical scenario (f3 = 1/8 and v = 1). 



Figure |5(b)| shows the scaling of the maximum of the 
Ising susceptibility Xa versus L. The slope 7/z/ = 1.76 ± 
0.02 is very close to the expected value of 7/4. We also 
plot the temperature T C {L) of th e the maximum of the 
specific heat at size L (Fig. |5(c) ), showing the expected 
asymptotic scaling T C (L) ~ + AjL x l v with v = 1. 
Finally, we computed the fourth-order Binder cumulant 



B(L,T) = 1- 



1 



3 la 



2\2 



(11) 



which shows a c harac teristic L-independent crossing at 
T c ~ 0.97 (Fig. |5(d)[) . Note however that the value of 



the cumulant at the crossing B(L = 72, T cross i ng ) ~ 0.65 
remains larger than the universal value of the 2D Ising 
model (0.6107). This discrepancy possibly originates 
from the fact that L is not significantly larger than the 
crossover length-scal e bey ond which the two peaks in 
P(E, T c ) merge (Sec. [II A ), making it uneasy to obtain a 
reliable estimate of the fourth moment of the distribution 
of chiralities. We also point out that similar discrepan- 
cies were observed in the simulation of other enwjpgent 
Ising systems with continuous degrees of freedomBET 
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FIG. 5: (Color online) Finite-size scaling at, or near, the 
transition temperature, a) Rescaled chirality (Eq. ^) ver- 
sus rescaled temperature for 40 < L < 80. b) Log-log plot of 
the maximum of the Ising susceptibility vs. L. c) Scaling of 
the maximum of the specific heat vs. 1/L (log-linear plot), d) 
Evolution of the Binder cumulant B(L,T) with temperature, 
from L = 16 to L = 88, from bottom to top. 



C. Z2 vortices in the Ising- 



model 



We now turn to the identification of Z2 vortices in the 
model defined in Eq. ||. The better-known case of an 
SO(3) ground state manifold (ferr omag netically frozen 
chiralities) was detailed in Section. II B Once the chiral 
degrees of freedom are relaxed, the ground state manifold 
is enlarged to 0(3) = Z2 x SO(3) and the definition of 
vorticity enclosed in lattice loops requires some caution. 
Let us consider two sites i and j and their associated el- 
ements of 0(3), Mi and Mj. As long as Mi and Mj have 
the same chirality, no extra difficulty arises. However, 
if they have opposite chiralities (determinants), then the 
mere existence of a continuous path connecting the two 
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elements in SO(3) is ill-posed, since they each belong to 
a different 5*0(3) sector of 0(3). Hence the computation 
of the circulation fl(£) makes sense only for loops C en- 
closed in a domain of uniform chirality. In particular, the 
computation of the Z2 vorticity on plaquettes located in 
the bulk of uniform domains follows the lines detailed in 



Section [IB without modification. 



The case of non-uniform plaquettes, sitting on a chi- 
ral domain wall, may be addressed in an indirect way. 
We consider a closed loop C that i) only visits sites with 
chirality cr, = +1, ii) encloses a domain of opposite chi- 
rality <7j = — 1. i) ensures that fl(C) is well defined. 
On the other hand one can always define unambiguously 
N V (C), the number of vortex cores on uniform plaque- 
ttes inside L. If the chirality were uniform inside L then 
il(£) = (— l)-^"^) would hold. However, when £ en- 
closes a domain with reversed chirality, an extra con- 
tribution arises on the right hand side, coming from the 
non-uniform plaquettes sitting on the domain wall, which 
are not accounted for by (— l)^^). Since fl(C) — ±1 we 
obtain Q(C) — e v ,(—l) Nv ^ where e w = ±1 acts as the 
topological charge of the chiral domain wall. As a result, 
the present workaround yields the total Z2 charge carried 
by the chiral domain wall, which is always well defined. 

shows a typical configuration near T c . 



Figure 7(a 



Again, vortices in the bulk of Ising domains are indicated 
by (red) bullets A thorough study of typical configura- 
tions reveals that most of these vortex cores are actually 
"paired" with a near by charged domain wall (not rep- 
resented in Fig. 7(a)). Once the charge of the walls is 
appropriately accounted for, the total measured charge 
is 1 indeed. Note that such vortex/charged- wall pairs 
feature two Z2 charges, but involve the creation of only 
one vortex core. Hence they are energetically favored 
compared to genuine pairs of Z2 vortices in the bulk of 
chiral domains. This is evidenced, for instance, by the 
proliferation of vortices at lower temperatures in the full 
Ising-lP 3 than in the RP 3 model (compare the upper 
panel of Fig.| with Fig.|6(b)|). 

To quantify the pairing effect of Z2 vortices with 
charged walls, we impose a domain wall by splitting in 
the system in two parts with frozen but opposite chi- 
ralit ies ( periodic boundary conditions are used). Fig- 
ure S(a) shows the excess density of vortices near the 
domain wall, compared to the density in the bulk of the 
domains, as a function of the distance to the interface 
for T = 1.3. In agreement with the trend observed in 
Fig. 7(a) , vortex/charged- wall pairs are clearly favored 
compared to vortex/vortex pairs in the bulk. Moreover, 
this effect is robust: the excess density of vortices is siz- 
able in a wide range of temperatures. Overall we antic- 
ipate that the same mechanism will prevail near more 
complex interfaces, such as those obtained in the full 
Ising-RP 3 model at equilibrium. 

For completeness, we computed the core vortex den- 
sity no (defined as the number of vortex cores on uniform 
plaquettes divided by the number of such plaquettes) as 
a function of temperature for different lattice sizes: the 
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FIG. 6: (Color online) a) Excess density of vortices vs. dis- 
tance to the domain wall at T = 1.3 (see text), b) Vortex 
density in the bulk of chiral domains vs. temperature T, from 
L — 16 to L = 88, from bottom to top. c) Maximum of O^p-J 
versus specific heat maximum maxT(C v ) for 16 < L < 88. 



increase in the v ortex density at the transition temper- 
ature (Fig. |6(b)| ) scales with the system size. This is 
apparent on the associated susceptibility dnn/dT whose 
maximum scale s with L like the maximum of the specific 
heat (Fig. 3(c) ). This is to be compared with the non- 



critical behavior of t he v ortex density in the MP 3 model 
discussed in Section 1IB where dnn/dT remains finite at 
all temperatures. 
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FIG. 7: (Color online) Thermalized configuration at T ~ T c 
of a system of size L — 72 (only a 45 x 45 square snippet 
is shown). Red bullets indicate elementary square plaquettes 
that i) host a I2 vortex core and ii) have all four Oi equal. 
Plaquettes with homogeneous up (resp. down) chiralities are 
represented by green (resp white) squares, a) the width of 
bonds is proportional to Bij (see text), b) Bonds are 

drawn only if Bij > 1/2. 



D. Spatial correlations of discrete and continuous 
fluctuations 



Figure 7(a) also reveals that the bonds standing across 
chiral domain walls {piOj = — 1) are also the most "dis- 
ordered", or thickest, ones. This strong correlation is all 
the more apparent in Figure |7(b)| , where only the highest- 



energy, most disordered, bonds, are represented, defined 
as > 1/2: they are rather scarce in the bulk of Ising 
domains. Hence the two kinds of fluctuations (associated 
to the discrete Ising variables and to the continuous 4D- 
vectors, respectively) are both localized close to Ising do- 
main walls. As a result, the energy barrier for the forma- 
tion of a domain wall is considerably lowered compared 
to the pure Ising model (with all 4D-vectors frozen in 
a ferromagnetic configuration) . This is evidenced by the 
low transition temperature T c = 0.97 for the Ising transi- 
tion observed in the present Ising-MP 3 model, compared 
to T c = 2.269 for the Ising model in 2D. 



E. Small modifications of the Ising-MP model and 
tuning of the nature of the phase transition 

In view of the peculiar finite-size, first-order-like, be- 
havior of the energy distribution shown in Figure |4J, we 
argue that the Ising-MP 3 model could be near a tricrit- 
ical point in some parameter space. This is consistent 
with the nature of the phase transitions observed in a 
number of related classical frustrated spin models with 
similar "contentjL-i.e. spin- waves, Z2-vortices and Ising- 
like chiralities Ma 

In support of our claim, we present two distortions of 
the original Hamiltonian @ which preserve the ground 
state symmetry, hence the nature of the excitations 
above, and show that they undergo a first-order phase 
transition. 

As it turns out, a simple change in the bond energy 



Eij = -<7i£7j(4({T i • Vj) 2 - a) 



(12) 



from a = 1 in the original Hamiltonian (||) to a = 1.75, 
is sufficient to drive the order-disorder transition of the 
Ising variables towards first order. This can be seen in 
Figu re H, where both the maximum of the specific hea t 



(Fig. 8(a) ) and that of the chiral susceptibility (Fig. 8(b) ) 
scale as ~ L 2 . Furthermore, contrary to the Ising-KP 3 
case, the energy probability distribution P(E 7 T = T c ) 
remains bimodal for all lattice sizes, with a minimum 
that gets more and more pronounced upon increasing the 
system size (not shown). 

Increasing a from a = 1 in (|lj) clearly decreases the 
energy gap for a chirality flip — ► — er.;. However, the 
associated modification of the entropy balance between 
the ordered and disordered phase is less obvious. 

Entropic effects are more explicit and better controlled 
in the following family of continuous spin models on 2D 
lattices: 



(13) 



Indeed, it was-shown that for large enough p (p 
for XY spinsc3 and p > 16 for Heisenberg spinscSE^I) 
these systems undergo a phase transition of the liquid- 
gas type. Obviously, tuning p does not change the energy 
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FIG. 9: (Color online) Energy distributions at the transition 
for model (111). 





(b) 



FIG. 8: (Color online) a) Scaling as L 2 of the maximum of the 
specific heat and b) of the chiral susceptibility with increasing 
system size L for a = 1.75 in Eq. 12l 



scale (Eij € [0, 1]), but increasing p gradually pushes the 
entropy toward the highest energies. 

Hence we tweak the Ising-KP 3 model in a similar way, 
with: 



Ei 



-c H a j {^v i -v j ) 2 ' p -1) 



(14) 



Once again, Monte Carlo simulations of the distorted 
model ( |l4l ) shows the signatures of a first-order transi- 
tion, for p as small as p — 2 (Fig. ||). 

Overall the previous two models illustrate our claim 
on the proximity of the Ising transition in the Ising-KLP 3 
model with a first-order transition. 



IV. ENTROPY OF OiOj BONDS 
MODEL ANALOGY 



POTTS 



Here we take another route and propose a simple qual- 
itative analogy to explain how the coupling of chiralities 
Oi to the continuous degrees of freedom Vi drives the Ising 
transition in the Ising-KP 3 model close to a first-order 
transition. 



FIG. 10: (Color online) Density of state g e (E) for a single 
bond in model (||). The two colors correspond to e = maj = 1 
(green) and e = — 1 (red). The lowest energy (E = —3 per 



bond) is attained when (vi 



1 and (TiCTw 



To that end we compute the density of states g(E) 
for a single bond It has two contributions g{E) = 

g + {E) + g~(E) depending on the value of e = UiUj-. 



9 e (E) 



d 3 u, 



d 3 Vj5 {E + e(4({Tj • v 3 ) 2 - 1)) 

d 3 ^ (£ + e(4(u°) 2 -1)) (15) 

where rotational invariance was used to fix Vj = 
[1,0,0,0]. Explicit evaluation of the integral gives 



s 3 Js 3 

,3 ; 

s 3 



5 e (-3 < eE < 1) = 



eE 



2tt V 1 - (E 



(16) 



Figure 1^ shows the evolution oig ± {E) with E. Interest- 
ingly, the density of state is remarkably small in the vicin- 
ity of the (ferromagnetic) ground state value (oiOj = 1, 
Vi ■ Vj = 1, and E = —3). However, if the two trihe- 
dra have opposite chiralities, the lowest energy state is 
attained for orthogonal vectors Hi ■ Vj = 0, with energy 
E = — 1, and the associated density of states diverges. 

This strong entropy unbalance between ordered (e = 1) 
and disordered bonds (e = —1) is,-similar to that of the 
well-known q— state Potts models.E3 In this model, each 
"spin" a can take q different colors, and the interaction 
energy is E = for neighboring sites {i, j) with the same 
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color Uj = o j , and E 
density of states 



1 otherwise, with the resulting 



5(0 < E < 1) = qd(E) + q{q - 1)6(E - 1). 

Hence "disordered" configurations indeed carry more 
weight than the ferromagnetic ground state. This en- 
tropy unbalance is seen to increase with q and in 2D it is 
known to eventually drive the order-disorder transition 
towards first order for q > 4rB 

A similar feature is obtained upon distorting the Ising- 
RP 3 model as in (|l4|). Indeed, the computation of the 
density of states g p (E) for p = 2 yields 



<£ =2 (-3 < eE < 1) 



1 y/2- A /l~7E 
27 (1-eS) 3 / 4 ' 



(17) 



which has the same features as in Fig. flC)| except that it di- 
verges as ~ a; -3 / 4 at E = ±1, instead of ~ x~ x l 2 . Hence, 
the distortion of the bond energy shown in ( |l4| ) essen- 
tially increases the entropy unbalance, which ulti matel y 
drives the Ising transition towards first-order (Sec. [II E] ) , 
much in the same way as in the q— states Potts model. 

To elaborate further on the proximity to a first or- 
der transition, it is useful to recall some results from the 
real-space renormalization group (RG) treatment o£-the 
g-states Potts model.Eil As noted by Nienhuis et in 
2D, conventional RG approaches give accurate results in 
the critical regime (q < 4) but inexplicably fail to pre- 
dict the crossover to a discontinuous transition for q > 4. 
The authors proposed that it originates from the usual 
coarse-graining procedure, by which a single Potts spin is 
assigned to a finite region in real space, using a majority 
rule. Intuitively, this brutal substitution becomes phys- 
ically questionable when there is no clear majority spin 
in the domain, a situation that is likely to occur when 
the number of colors q is large enough. In particular, it 
yields the possibility of a ferromagnetic effective inter- 
action between such (artificially) polarized super-cells, 
even if the microscopic spins are disordered. Hence, con- 
ventional coarse-graining overestimates the tendency to 
ferromagnetic order. 

In Ref E it is argued that disordered regions interact 
only weakly with their neighbors, hence they are better 
coarse-grained as a vacancy (missing Pott spin). In sup- 
port to this intuitive picture, it was shown that, once 
the parameter space of the original Potts Hamiltonian is 
enlarged to include the fugacity of these vacancies, the 
real-space RG treatment of the Potts model is able to de- 
tect the crossover of the order-disorder transition, seen as 
a liquid-gaz transition of the vacancies. 

In the following section we discuss a simplified version 
of the Ising-MP 3 model where discrete variables tj are 
introduced. The latter play a role similar to that of the 
vacancies in Ref. BlL 



V. EFFECTIVE DILUTED ISING MODEL 

In this section we propose a simplified model, with 
strong analogies to the Ising-IRP 3 model (||), that cap- 
tures the spatial correlations evidenced in Fig. |7(b)| and 
where the entropy unbalance discussed above for the 
Potts model, is at play. 

We replace the vector degrees of freedom with discrete 
variables tj = 0, 1, so that the energy becomes: 

E = -J2 Wi&Wj - !) + D ( T ) U - ( 18 ) 

We introduce a temperature-dependent "chemical poten- 
tial" D(T) to tune (U) (hence (titj)). The relation with 
the original model (g) can be understood as follows. A 
site with U = 1 represents a vector which is collinear (or 
almost collinear) with the "majority" of its neighbors. 
On the other hand, a site with tj — represents a vec- 
tor which is perpendicular (or almost perpendicular) to 
the majority of its neighbors. The fact that the vector- 
vector correlation length is significantly larger than one 
lattice spacing at the temperatures of interest justifies 
that, locally, the vectors have a well-defined local ori- 
entation. Then, we simply replace (?7j ■ Vj) 2 by titj. Of 
course, in the original model, two vectors Vi and Vj can be 
simultaneously i) orthogonal to most of their neighbors 
(fi = tj = 0) and ii) parallel to each other ((vi-Vj) 2 = 1): 
such situations are clearly discarded by this discretized 
model.c3 As a final encouragement to study the discrete 
model of Eq. |l8|), we mention that its single bond den- 
sity of state is qualitatively similar to that of the original 
model (Fig. |l^): it has two ground states at E = — 3 
(uj = Uj and ti = tj = 1), and 2x3 excited states asso- 
ciated to chiral flip at E = — 1 (<7j = — Oj and titj = 0). 

This modei-jclosely resembles the celebrated Blume- 
Capel modelj 3 ^! where an Ising transition becomes first 
order when the concentration of "holes" (sites with tj = 0) 
becomes large enoughO Since a simple mean-field ap- 
proximation is sufficient predict the first and second or- 
der transition lines of the Blume-Capel model has (de- 
pending on the crystal field parameter A) we deter- 
mine the mean-field phase diagram of Eq. [18|, using two 
mean-field parameters (ti) = t and (crj) = a. Fig. |ll| 
shows that it is composed of four transition lines, and 
one obtains four distinct order-to-disorder transitions 
depending on the value of D(T). Namely, upon in- 
creasing D(T) we find i) a second order ferromagnetic- 
paramagnetic transition (large negative D) , ii) a first or- 
der ferromagnetic-paramagnetic transition, iii) a first or- 
der ferromagnetic-antiferromagnetic transition and iv) a 
second order antiferromagnetic-paramagnetic transition. 
Hence we conjecture that the present model also has a tri- 
critical point, where the ferromagnetic to paramagnetic 
transition changes from second order to first order. This 
approach is obviously too crude to be quantitatively ac- 
curate, but we can still make contact with the Ising-KLP 3 
model by adjusting the chemical potential D(T) to en- 
force (Utj) = ((vi ■ Vj) 2 ) (the right-hand side is com- 
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D 



FIG. 11: (Color online) Schematic mean-field phase diagram 
of the discrete model ( |18| ) in the (T, D) plane. The two thin 
lines are trajectories of our models (|l4), with p = 1,2). 



Formally, the integration over the 4D-vectors leads to 
the following energy E e s for a configuration {cr,} of the 
chiralities: 

£ cff (W) = -Tin (( e -^«"«.*i»)) , (19) 

where E({rji,Vi}) is given by Eq. |[ and (•••) is the 
T = oo average for each vector spin Hi € S 3 (uniform 
measure on S 3 x • ■ ■ x S 3 ). 

In the following, we derive the effective interaction of 
the chiralities by expanding Eq ( |l9| ) in powers of (3 = 1/T 
up to order f3 a . 

B. 1/T expansion 



puted in the original Ising-MP 3 model). Upon chang- 
ing the temperature, this model describes a curve in the 
D — T plane such as those shown in Fig. [ll]. At T = oo, 



and one can show that D(T = oo) 



0. 



Further, upon decreasing the temperature ((«, ■ Vj) ) 
decreases and D(T) decreases, until the ferromagnetic 
phase of the discrete model is reached. 

Overall this supports our claim that the unusual finite- 
size behavior of the Ising-MP 3 model (Fig. |j), as well 
as t he pro ximity to a first order chiral transition (Sec- 
tion [HE) originates from a nearby tricritical point in 
parameter space. Further, the present study suggests 
that the nature and origin of the tricritical point can 
be understood, at least qualitatively, from an effective 
Blume-Capel like model. 



VI. EFFECTIVE ISING MODEL WITH 
MULTIPLE-SPIN INTERACTIONS 

In this section we detail a more quantitative approach 
to the Ising-MP 3 model, in which the vector spins Vi are 
integrated out perturbatively in = 1/T, in order to 
derive an effective Ising model for the chirality degrees 
of freedom. From this standpoint, multiple spin interac- 
tions between chiralities are directly responsible for the 
"proximity" to a first order transition. 



A. Integrating out the vector spins 

Because the order parameter of the transition is the 
chirality, and because the vector spins never order at 
T > 0, it is natural to look for an effective model involv- 
ing only the chiralities. Moreover, we have shown that 
the Ising domain walls are accompanied by short distance 
rearrangements of the 4D-vectors, and that it is a very 
important aspect of the energetics of the system. It is 
thus natural to expect that a high temperature expan- 
sion for the vector spins (which captures short-distance 
correlations) will be semi-quantitatively valid. 



x E E ••• E 



y> (-/?)" 
(Ei 1 j 1 Ei 2 j 2 ■ ■ ■ Ei n:j 



with 



Ej 



l-4(*.tf,) ! 



(20) 



(21) 



As (Eij) = 0, expanding the logarithm of Eq. |9|in pow- 
ers of j3 yields the following cumulant expansion: 



/? 3 



of E E C iijii 2 h a ii a h^h 



2! 



(»l,Jl> (*2,J2> 



3! 



E E E 



(7 il (T jl o- i2 a h o- i3 (Tj 



(«2J2> (l3,j3) 



+ 4! E "" E ^iji-i4i4 f7 «i f7 ii 

ih-jl) (l4j4> 



e>(/3 5 



with the cumulants 



CI 



(E lljl E i . 



C 



»1J1^2J2*3J3 



(Ei 1 j 1 Ei 2 j 2 Ei 3 j 3 ) 



(22) 

(23) 
(24) 



~ (Ei 1 j 1 Ei 2 j 2 ) (Ei 3 j 3 Ei 4 j 4 ) 

— (Ei 1 j 1 Ei 3 j 3 ) (Ei 2 j 2 Ei 4 j 4 ) 

- (E njl E iij4 ) (E i2j2 E l3j3 ) (25) 



As usual in series expansion, the cumulant C" is 
non-zero only if the graph defined by the n bonds 
• • • ,(in,jn) is connected. Moreover, rotational 
invariance ensures that only graphs that are one-particle- 
irreducible contribute. For this reason, C? - ,• ,• ls non ~ 

< 1J1 '2J2 

zero only when the two bonds coincide, resulting in a 
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constant contribution (independent of the Oi) to E c g. At 
the next order, and for the same reason, the only non- 
zero C come from graphs where the three bonds coin- 
cide (and C121212 = !)• This generates an effective first 
neighbor Ising interaction proportional to 1/T 2 : 

C 4 only provides a constant contribution to the effec- 
tive energy. C 5 and C 6 terms introduce new two-chirality 
interactions, between first, second and third neighbors. 
Moreover, a C 6 term gives the first interaction with more 
than two chiralities, namely: 

"^5 Yj g ^3 g ^u (27) 

(i,j,k,l) 

where the sum runs over square plaquettes. 

At this order in f3, the Ising-JRi-3 model appears as an 
Ising model with two and four spin interactions. The ef- 
fect of such multiple-spin interactions has been studied 
for the 3D Ising model, where an additional four spin in- 
teraction, if_.it is large enough, can make the transition 
first order£3 This can be understood, at least qualita- 
tively, from a very simple mean field point of view. In- 
deed, p— spin interactions will translate into terms of the 
order of mP in the Landau free energy (to being the order 
parameter) . Hence it is clear that tuning the strength of 
multiple (> 4) spin interactions can reshape the free en- 
ergy landscape, and drive the transition from second to 
first-order. 

However, various approaches predicted that the sim- 
plest 4-spin interactions were not enoueh-,te obtain a 
first-order transition in two dimensions BJcMa To check 
these predictions we performed Monte Carlo simulations 
of a simple Ising model with first neighbor coupling, sup- 
plemented with a 4-spins plaquette interactions The 
Hamiltonian reads 

H = - J ^ a i a 3 ~ A'^CrjCTjCTfcCTi, (28) 

<ij> □ 

with J,K > 0. Using finite-size scaling analysis, we find 
that the transition remains of second-order on the whole 
range < K/J < 10. 

This lead us to continue the high-temperature expan- 
sion up to order /3 8 , where a multi-spin interaction in- 
volving 6 chiralities is generated (with new 2- and 4- 
interaction terms). The effective Hamiltonian becomes 
quite complicated and therefore we resort to a simple 
mean-field calculation. 



C. Mean-field approximation 

The expansion of Eq. [ll] up to /3 8 leads to a huge num- 
ber of terms and it is necessary to proceed systemati- 
cally in order to obtain all the diagrams. In mean-field, 



TABLE I: List of all diagrams contributing to the high tem- 
perature series expansion up to order /3 8 . All vertices with an 
odd number of lines correspond to a Oi in the effective Ising 
model, to a m in the mean-field approximation. M g is the 
number of ways in the graph g can be located on the lattice. 



each chirality is replaced by its mean value (en) = to, 
which simplifies the diagrammatic expansion of the ef- 
fective Hamiltonian. We have written a (Maple) sym- 
bolic code that i) generates all possible diagrams on the 
square lattice, ii) assigns a weight to to multiply con- 
nected vertices with an odd number of bonds iii) com- 
putes the integrals over the different vectors exactly and 
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iv) computes the number of ways J\f g the graph g can 
be located on the lattice. Af g corresponds to the prod- 
uct of the number of bond ordering by the number of 
transformations (rotations, reflections, deformations on 
the lattice) changing the representation of the graph but 
not the graph itself. It is also worth noting that graphs 
free of articulation vertex have zero cumulants. In ad- 
dition, we discard graphs that give irrelevant constant 
contributions, i.e. graphs where the multiplicity of each 
vertex is even. With these prescriptions, only 55 graphs 
contribute at order (3 8 . They are listed in Tab|j. 

As a result, we obtain the effective energy in the mean- 
field approximation: 



E eff (m) 



(3 2 7/3 4 
3 45 
5173/3 8 \ 2 
145800 J 
(3 5 7/3 7 
18 ~ 162 

81 



(3 5 
~9~ 



391/3 6 
4536 



81 



2/?% 4 



(29) 



Hence the mean-field free energy is given by 
F(m) = E(m) — TS(m) where the entropy 
S(m) = -i±!2 1n(±4p) - i^ln(i=^). F(m) is 
minimized with respect to the magnetization m: a phase 
transition occurs at T = 1.07 between an ordered phase 
m^O and a paramagnetic phase m = 0. Further, the 
transition is of first order albeit with a very small free 
energy barrier (see Fig. |l2|) . 

This computation is repeated for the Ising-Rp^ model 
with "non- linear" interaction, (Eq. |lj with p = 2). We 
obtain 



f3 2 81/? 4 1 3/3 5 1281493/3 6 



E eS (m) = 1 - 



3 800 200 
4877/3 7 5746857169/? 6 
80000 ~ 

13/3 5 



10080000 



+ - 



82944000000 
1351/3 7 227/3 8 



400 



32000 



8000 



227/3 8 \ 
' 48000 J 



(30) 



As can be seen in Figure [T^, the first-order transition in 
this model is much stronger than in the p = 1 case, i.e. 
(A/3 c P) p=2 3> (A/3 c F) p =i. We also mention that similar 
conclusions can be made for model jl^ ) with a = 1.75. 

Once again, these results support our claim that the 
Ising-KP 3 model is close to a point in parameter space 
where the transition changes from first to second order, 
in qualitative agreement with the simulation results. 
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FIG. 12: (Color online) Mean-field free energy for model ( p^ ) 
at the transition temperature, for p = 1 (Ising-RP 3 model, 
solid line) and p = 2 (dashed line). In both cases there are two 
stable mean-field solutions and the transition is first order. 
However the energy barrier is much smaller for p = 1. 



VII. CONCLUSION 

We have proposed a minimal Ising-RP 3 model that 
captures most of the low energy features of frustrated 
Heisenberg models with an 0(3)— manifold of ground 
states. The discrete symmetry-breaking associated with 
the existence of two connected components of 0(3) leads 
to an Ising-type continuous transition with unconven- 
tional features at small sizes. The complexity of this 
transition is due to the strong coupling of the Ising chi- 
rality fluctuations to short range continuous spin fluctu- 
ations and to the topological defects of these textures. 

We provided a consistent picture of these Z2 defects 
in presence of chiral fluctuations. Namely, we showed 
that chiral domain walls can also carry a Z2 topolog- 
ical charge, albeit delocalized over the interface. In- 
spection of configurations revealed that isolated vortex 
cores are almost always located nearby charged domain 
walls. At the transition temperature, most defects con- 
sist in vortex/charged- wall pairs, while there are a few 
vortex/ vortex pairs and almost no isolated defect. This 
mechanism is essential in understanding the appearance 
of point defects at a temperature much lower than in the 
case where chiral fluctuations are absent. 

We also studied some variants of the model and showed 
that the continuous transition easily becomes first-order, 
which lead us to conjecture the existence of a nearby tri- 
critical point in parameter space. We clarified the role 
of short-range fluctuations of the continuous variables 
in this mechanism by an analogy with the large g-state 
Potts model. In this analogy the high density of states 
with orthogonal 4D-vectors translate into the large num- 
ber q(q — 1) of states with disordered bonds of the Potts 
model. Finally we studied two derived versions of the 
original model: i) a diluted-Ising model, which somehow 
corresponds to a coarse-grained version of the Ising-KLP 3 
model (the continuous variables are locally averaged and 
replaced by discrete variables), and ii) an effective mul- 
tispin Ising model, obtained by tracing out the vector 
spins, order by order in (3 (up to (3 8 ). These two models 
predict either a weakly first-order or a continuous tran- 



13 



sition, as well as the existence of a tricritical point. 

This study sheds light on the large variety of behav- 
iors reported rfeftrfhiral phase transitions in frustrated 
spin sy stems. QlHjI3 In all these models, the chirality is 
an emergent variable more or less coupled to short-range 
spin fluctuations: in the J\ — J3 model on the square 
lattice the chiral variable is the pitch of an helix, the 
coupling to the short range spin fluctuations is probably 
small and the transition appears to be dearly second- 
order and in the Ising universality class In the cyclic 
4-spin exchange model on the triangular lattice, the or- 
der parameter is a tetrahedron and the chiral variable 
is associated to the triple product of three of these four 
spina,— .the transition is probably very weakly first or- 
der.aE3 In the J\ — J2 model on the kagome lattice, the 
order parameter at T = is a cuboactedron, and the 
chiral variable is associated to the triple product of three 
of these twelve spins. The phase transition evolves from 
weakly to strongly first order when tuning the parameters 
towards a ferromagnetic-antiferromagnetic phase bound- 
ary at T = 0: this can be understood in the light of 
the present work. Tuning the parameters towards the 
ferromagnetic phase frustrates the 12-sublattice Neel or- 
der and favors short-range disorder and vortex formation. 
The associated increase in the entropy unbalance drives 
the transition from Ising to strongly first order, much in 



the same way as in Section IV. 

On the technical side, this proximity of a tricritical 
point in parameter space evidences why simulations and 
experiences must be lead with great caution, a conclusion 
equally supported by recent work from the quite different 
standpoint-of the non perturbative renormalisation group 
approach. Ell 
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APPENDIX A: MONTE-CARLO ALGORITHM 

In rtius section we detail the Wang-Landau algo- 
rithmE3'E3 used to simulate model (Q). This method 
consists in building the density of states g(E) progres- 
sively, using successive Monte Carlo iterations. Elemen- 
tary moves consist in rotations of vectors as well as 
flips of chiralities <7j. In the case where the chiralities 



are frozen (Section II B ) only the rotation movements 
are performed. For completeness we mention that the 
four-vectors V[ are sampled uniformly on S 3 using 3 ran- 
dom numbers (r, 77 and u), independent and uniformly 



distributed in [0, 1], according to 

fr cos(2nr)) 
fr sin(27r?7) 



y/1 — r cos(2-7w) 
_ yl — r sin(27w) 



(Al) 



Starting from an initial guess g(E), the acceptance of a 
trial flip/rotation is decided by a Metropolis rule 



II(o -> n) = Min 1, 



9{Eq) 

g(E n ) 



(A2) 



where the subscripts o and n corresponds the old and 
new configurations, respectively. 

Every time a configuration with energy E is visited, 
the density of states g(E) is multiplied by a factor / > 1, 
g(E) <— g(E)f. To ensure that all configurations are 
well sampled, a histogram H(E) accumulates all vis- 
ited states. The first part of the run is stopped when 
H(E) > 10 4 . In a second part, the histogram H(E) is 
reset and the run is continued, but the modification fac- 
tor / is-jaow decreased to fi < /. In the original paper by 
Wang,E3 fi was taken as \ff, but this .choice is not nec- 
essarily the best for continuous models .a In this case the 
convergence properties were found to be quite satisfying 
upon choosing /1 = f 07 . The random walk is contin- 
ued until the histogram of visits H(E) has become "flat". 
Once again, H(E) is reset and the modification factor is 
decreased to < /1, etc... Accurate density of states 
g n {E) are generally obtained at the n-th iteration, where 
n is such that /„ is almost 1 (typically /„ — 1 < 10 -8 ). 

Since this model has continuous variables, its energy 
spectrum is also continuous, and the choice of the energy 
bin requires special care: if the bin is too large, important 
details of the spectrum may be lost, whereas if it is too 
small, a lot of computer time is wasted to ensure the 
convergence of the method. In the temperature range 
of interest, a size of order ~ 0.1 is a good compromise 
for the model of Eq. El. In addition, the energy range 
is limited to the region relevant at the transition, and 
yields all thermodynamic quantities accurately at these 
temperatures. 

Once the density of states g(E) is obtained, all mo- 
ments of the energy distribution can be computed in a 
straightforward way as 



(E n ) = 



fg(E)E n exp(-[3E) 



(A3) 



fg(E)exp(-PE) 
from which the specific heat per site is readily obtained: 



Cv = „ ((E 2 ) - ((E)) 2 ) 



(A4) 



For the thermodynamic quantities that are not directly 
related to the moments of the energy distribution, such 
as the chirality, the vorticity, and their associated sus- 
ceptibilities, or the Binder cumulants, we proceed as fol- 
lows: an additional simulation is performed where g(E) 
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is no longer modified (a "perfect" random walk in energy 
space if the density of states g(E) is very accurate). In 
this last run additional histograms are stored: chirality 
histograms cr(E), a 2 (E), a 4 (E), and vorticity histogram 
V(E), V 2 (E), V\E). 

Chirality (or vorticity) moments are then obtained 
from the simple ID integration 



(AT 



fg(E)exp(-(1E) 



(A5) 



from which, say, the Binder cumulant, is readily obtained 



as 



U = 1- 



(M 4 ) 
3((M 2 )) 5 



(A6) 



where M = o or ^—.Contrary to the method of the Joint 
Density of States Jlil our method does not require the 
construction of the two dimensional histogram g{E,cr). 
Hence it is not limited to modest lattice sizes (the above 
quantities are computed for L up to L = 80). 



J. Villain, J. Phys. (Paris) 38, 385 (1977). 
P . Chandra. P. Coleman, and A. I. Larkin, 



tt. 64, 88 (1990). 



s. Rev 



25 
26 



Rev. Mod. Phys. 51, 591 (1979' 



N. D. Mermin, 
T. Mom oi, H. Sakamoto, K. Kubo, |Phys. Rev. B 59, 9491 



(1999) 



C. Weber. L. Capriotti, G. Misguich, F. Becca, M Elhajal 



aid F. Mila, 



Phys. Rev. Lett. 91, 177202 (2003) 



E. Domany, M. Schick , and R. H. Swendsen, Phys. Rev 



T. Momoi, K. Kubo, and K. Niki, [Phys. Rev. Lett. 79 



2)81 (1997) 



Lett. 52, 1535 (1984) 



H. W.J. Blote, W. Guo, and H. J. Hilhorst, 



Lett. 88, 047203 (2002) 



Phys. Rev 



L Capriotti and S. Sachdev, Phys. Rev. Lett. 93, 2572061 ZM A. C. D. van Ente r and S. B. Shlosman, [Phys. Rev. Lett 



( ^004) 

J. -C. Domenge, P. Sindzingre, C . Lhuillier and L. Pierre, 



Phys. Rev. B 72, 024433 (2005) 



31 



89, 285702 (2002) 



J. Phys. C 6, L 445 (1973) 



J. Baxter, 

Nienhuis. A. N. Berker. E. K . Riedel and M. Schick 



J.-C. D omenge, C. Lhuillier, L. Messio , L. Pierre and 



PViot, Phys. Rev. B 77, 172413 (2008) 



Phys. Rev. Lett. 43, 737 (1979) 



P. C. Hohenberg. Ph ys. Rev. 158. 383 (1967): N. D. Me r 



M. Blume, Phys. Rev. 141, 517 (1966); H. W. Capel, 



min and H. Wagner 

s. 



Phys. Rev. Lett. 
joyeul Lee and Koo-ChrTTee~ "]p*h'vs 




17. 1133 (1966) 



s. Rev. B 49. 15184 



Physica 32 966, (1966) 
O. G. Mouritsen, S. J . Knak Jensen, B. Frank, [Phys 



Rev. B 24 



D. Loison and P. Simon, Phys. Rev. B 61, 6114 



E. Bre zin and J. Zinn- Justin, Phys. Rev. Lett. 36, 691 



C.976) 



Mukamel 
J. Oitmaa and R 



347 (1981): O. G. Mouritsen. B . Frank and D. 



Phys. Rev. B 27, 3018 (1983) 



C: Solid State 



W. Maier and A. Saupe, Z. N aturforsh. A 14, 882 (1959) . 
P. A. Lebwohl and G. Lasher, |Phys. Rev. A 6, 426 (\V72\ 
P . E. Lammert, D. S. Rokhsar and J. Toner, [Phys. Rev 



W. Gibberd, J. Phys 
Phys. 6, 2077 (1973). 

J. Oitmaa, J. Phys. C: Solid State Phys. 7, 389 (1974). 
M. Gitterman and M. Mikulinski, J. Phys. C: Solid State 
Phys. 10, 4073 (1977). 

B. Delamotte, D. Mouhanna and M. Tissier, 



tt. 70, 1650 (1993). P. E. Lammert, D. S. Rokhsar and 



J Toner, phys. Rev. E 52, 1778 (19951 
S Car acciolo. R. G. Edwards. A. Pelis setto and A. D. 



B 69, 134413 (2004) 



Phys. Rev 



F. W. Wang and D. P. Landau, Phys. Rev. Lett. 86, 2050 



Sokal, |Phys. Rev. Lett. 71, 3906 (1993j 



(2001) 



H. Kawamu ra and S. Miyashita, J. Phys. Soc. Jpn. 53 



4 .38 (1984) 



B. J. Schulz. K. Binder. M. Muller and D. P. Landau, 



Phys. Rev. E 67, 067102 (2003) 



P. Poulain . F. Calvo. R. Antoine. M. Br over and Ph 



G. Kohring and R. E. Shrock, Nucl. Phys. B 285, 
(1987). 



504 



Dugourd, [Phys. Rev. E 73, 056704 (2006) 



41 C. Zhou. T. C. Schulthess. S. Torbrug ge and D. P. Landau 



T. Dombre and N. Read, [Phys. Rev. B 39, 6797 (1989' 



18 
19 



Azaria. B. Delamo tte. and T. Jolicoeur, Phys. Rev 



P. 
L 

H. Kunz and G. Zumbach. 
P. Azaria, B. Delamotte, T 



tt. 64, 3175 (1990) 



Phys. Rev. Lett. 96, 120201 (2006) 



Phys. Rev. B 46, 662 (1992). 



Joli coeur, and D. Mouhanna, 



Phys. Rev. B 45, 12612 (1992) 



B. W. Southern and A. P. Young, Phys. Rev. B 48, 13170 



993) 



B W. Southern and H.-J. Xu, Phys. Rev. B 52, R3836 



( 1995) 



M. Hasenbusch, |Phys. Rev. D 53, 3445 (1996)|. 

M. Wintel, H . U. Everts and W. Apel, |Phys. Rev. B 52 
1 ^480 (19951 

M. Caffarel, P. Azaria, B. Delamotte, and D. Mouhanna 



Phys. Rev. B 64, 014412 (2001) 



J. Phys. A 9, 2125 (1976). 

42 This model does not encompass the breathing modes of 
the trihedra that are present in the full frustrated spin 
systems: we suspect that these modes would not change 
qualitatively the present picture. 

43 jg>p« j s £ ne rea j p ro j ec tive space. It is formed by taking the 
quotient of R" +1 — {0} under the relation of equivalence 
x ~ Xx for all real numbers A 7^ 0, or equivalently, by iden- 
tifying antipodal points of the unit sphere, S n , innR. n+1 . 
RP 2 was introduqeA4a the context of liquid crystalst-H (see 
also Refs. 0,p.laJ 

44 Overall our simulations show no obvious signature of a 
vortex unbinding transition. However, more conclusive ar- 
guments require the computation of more vortex-sensitive 
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quantities, such as then spin stiffnesses, or area versus 
perimeter scaling lawsJI-J which are numerically very de- 
manding. 

In the Ising-RP 3 model, domain walls separate regions 
with a — + 1 from those with a = — 1. Fig. |icj shows that 
ferromagnetic configurations of the vectors are favored in 
the bulk of chirality domains. On the other hand, orthog- 
onal configurations are favored on bonds that stand across 
domain walls. Hence the lowest energy configuration for 
the vectors is to point in direction, say, [1,0,0,0] in the 
a = +1 domain, and in some orthogonal direction, say 
[0,1,0,0] in the a — — 1 domain. An analog low-energy 
configuration in the discrete model of Eq. [ISfl is obtained 



by setting U = 1 in the bulk of both domains, and tj = 
in the vicinity of the domain wall. i — i 
The Hamiltonian of the Blume Capel model 3 -^ is H = 
-Jj] <j j > SiSj + AJ] 4 S?, where Si denotes a spin-1 at 
site i, J the ferromagnetic interaction and A the crystal- 
field. The Blume-Capel model and the site-diluted spin 
model @ become equivalent once the —1 contribution is 
dropped in Eq. [l8| and D(T) = A - Tln(2). 
Although the two-spin interactions appearing in E e g up to 
order T~ 6 are not strictly limited to first neighbors (see 
the graphs in Table ||). 



